;;......................................................
;; Prepare fire emission data for CAM5
;;
;; Input files:
;;    GFED daily emission data
;;    CAM5 elevated OC data (for profile)
;;
;; Author : Kai Zhang (kai.zhang@pnnl.gov)
;;......................................................

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin

   yy = getenv("yy")
   tp = getenv("tp")

   ;;...................................................... 
   ;; Input files 
   ;;...................................................... 
   fna = "GFED3.1_daily_"+tp+"_"+yy+"_2x2.5.nc" 

   if(tp.eq."BC") then 
      fnb = "ar5_mam3_bc_elev_2000_c090726.nc" 
   end if 

   if(tp.eq."OC") then 
      fnb = "ar5_mam3_oc_elev_2000_c090726.nc" 
   end if 

   if(tp.eq."SO2") then 
      fnb = "ar5_mam3_so2_elev_2000_c090726.nc"
   end if 

   fnc = "ar5_mam3_num_a1_elev_2000_c090726.nc"

   ;;...................................................... 
   ;; set date 
   ;;...................................................... 
   load "date_2008.ncl" 

   if(yy.eq."2005") then 
      date = date - 30000
   end if 

   if(yy.eq."2006") then 
      date = date - 20000
   end if 

   if(yy.eq."2007") then 
      date = date - 10000
   end if 

   if(yy.eq."2009") then 
      date = date + 10000
   end if 

   if(yy.eq."2010") then 
      date = date + 20000
   end if 

   if(yy.eq."2011") then 
      date = date + 30000
   end if 

   ;;printVarSummary(date) 
   print(" ") 
   print(" ") 

   fla = addfile(fna,"r") 
   flb = addfile(fnb,"r") 

   ;;...................................................... 
   ;; emission date 
   ;;...................................................... 
   print("read data from :    "+fna) 
   print(" ") 
   print(" ") 

   va = fla->$tp$
   tima = fla->time

   lon = flb->lon
   lat = flb->lat

   ;;...................................................... 
   ;; cam5 data 
   ;;...................................................... 
   print("read data from :    "+fnb) 
   print(" ") 
   print(" ") 

   altitude = flb->altitude
   altitude_int = flb->altitude_int

   nt = dimsizes(tima) 
   nlev = dimsizes(altitude) 
   nlat = dimsizes(lat) 
   nlon = dimsizes(lon) 

   ;;...................................................... 
   ;; total fire emission of cam5
   ;;...................................................... 
   vv = flb->forestfire ;;(time, altitude, lat, lon)  
   vu = flb->grassfire ;;(time, altitude, lat, lon)  

   vv = vv + vu 

   ;;...................................................... 
   ;; profile 
   ;;...................................................... 
   debug = False ;; True

   if(debug) then 
     print(" Jan Emission Profile before remapping ") 
     print(" ") 
     print(" "+vv(0,:,70,50)) 
     print(" ") 
     print(" ") 
     print(" Nov Emission Profile before remapping ") 
     print(" ") 
     print(" "+vv(10,:,70,50)) 
     print(" ") 
     print(" ") 
   end if

   vsum = dim_sum(vv(time|:,lat|:,lon|:,altitude|:)) 

   do ilev = 0,nlev-1    
      vv(:,ilev,:,:) = vv(:,ilev,:,:) / (1.e-33 + vsum(:,:,:) ) 
   end do 

   if(debug) then 
     print("====================") 
     print(" Column Total Emission Flux ") 
     print(" ") 
     print(" Jan: "+vsum(0,70,50)) 
     print(" Nov: "+vsum(10,70,50)) 
     print(" ") 
     print(" ") 
     print("====================") 
     print(" Jan Ratio Profile before remapping ") 
     print(" ") 
     print(" "+vv(0,:,70,50)) 
     print(" ") 
     print(" ") 
     print("====================") 
     print(" Nov Ratio Profile before remapping ") 
     print(" ") 
     print(" "+vv(10,:,70,50)) 
     exit
   end if

   ;;...................................................... 
   ;; time
   ;;...................................................... 

   time = new(nt,"double") 
   time = ispan(1,367,1) - 2 
   time@calendar = "Gregorian"
   time@units = "days since "+yy+"-01-01 00:00:00" 
   time = ut_convert( time, "days since 2005-01-01 00:00" ) 
   time@long_name = "Time"

   ;;...................................................... 
   ;; date 
   ;;...................................................... 

   date!0 = "time" 
   date&time = time 
   date@units = "YYYYMMDD"
   date@long_name = "Date"

   ;;...................................................... 
   ;; output data 
   ;;...................................................... 

   print("define new variable :      " ) 
   print(" ") 
   print(" ") 
 
   vb = new((/nt,nlev,nlat,nlon/),"float") 

   vb!0 = "time" 
   vb&time = time
   vb!1 = "altitude" 
   vb&altitude = altitude
   vb!2 = "lat" 
   vb&lat = lat
   vb!3 = "lon" 
   vb&lon = lon

   vb@long_name = tp+" fire emissions from GFED v3.1"
   vb@units = "molecules/cm3/s" 
   vb@_FillValue = -999.

   ;;...................................................... 
   ;; vertical profile remapping 
   ;;...................................................... 

   idx = (/1,  32, 60, 91,  121, 152, 182, 213, 244, 274, 305, 335 /) 
   idy = (/31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 365 /) 

   idd = idy - idx 

   print("nday in mon : "+ idd) 
   print(" ") 
   print(" ") 

   do im = 0, 11 
      ida = idx(im)  
      idb = idy(im) 
      print("imon, ida, idb, id : "+im+" "+ida+"   "+idb) 
      vb(ida:idb,:,:,:) = conform(vb(ida:idb,:,:,:), vv(im,:,:,:), (/1,2,3/) ) 
   end do

   vb(0,:,:,:)   = (/ vv(11,:,:,:)  /) 
   vb(366,:,:,:) = (/ vv(0,:,:,:)  /) 

   do ilev = 0,nlev-1   
      vb(:,ilev,:,:) = (/ vb(:,ilev,:,:) * va(:,:,:) /) 
   end do 

   ;;...................................................... 
   ;; emission layer height (cm) 
   ;;...................................................... 

   dalt = altitude 

   do ik = 0, 13-1
      dalt(ik) = ( altitude_int(ik+1) - altitude_int(ik) ) * 1000. * 100. 
   end do 

   print(" ") 
   print(" ") 
   print("Emission Layer Height : "+dalt) 
   print(" ") 
   print(" ") 

   vc = vb 

   vc(:,:,:,:) = conform(vc(:,:,:,:), dalt(:), (/1/) )

   vb = vb / vc 

   ;;...................................................... 
   ;; unit (g/m2/day -> molecules/cm2/s) 
   ;;...................................................... 

   mole_to_molecules = 6.022e23

   m2_to_cm2 = 1.e4 

   day_to_s = 86400. 

   vb = vb * mole_to_molecules / m2_to_cm2 / day_to_s 

   if(tp.eq."BC") then 
      vb = vb / 12. 
   end if 

   if(tp.eq."OC") then 
      vb = vb / 12.
   end if 

   if(tp.eq."SO2") then 
      vb = vb / 64.
   end if 

   ;;...................................................... 
   ;; output data 
   ;;...................................................... 

   fno = "emis_GFED3.1_daily_"+tp+"_"+yy+"_2x2.5.nc" 

   print("output file : "+fno) 
   print(" ") 
   print(" ") 
 
   system("rm "+fno) 

   flo = addfile(fno,"c") 

   varname = tp+"_fire"

   flo->$varname$=vb 
   flo->lat=lat
   flo->lon=lon
   flo->date=date
   flo->altitude=altitude
   flo->altitude_int=altitude_int

   ;;...................................................... 
   ;; create global attributes of the file
   ;;...................................................... 
   fAtt               = True          
   fAtt@Author        = "Kai Zhang (kai.zhang@pnnl.gov)"
   fAtt@Source        = "GFED3.1 fire emission data and CAM5 input files"
   fAtt@creation_date = systemfunc ("date")
   fileattdef( flo, fAtt )

end


